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| SUMMARY 

I 

i 

Methods are described by which the dynamic interaction of structures 
With surrounding fluids can be computed by using finite element techniques. 

In all cases, the fluid is assumed to behave as an acoustic medium and is 
initially stationary. Such problems can be solved either by explicitly 
modeling the fluid (using pressure or displacement as the basic fluid unknown) 
jor by using decoupling approximations which take account of the fluid effects 
without actually modeling the fluid, 
f 

! INTRODUCTION 


Recently there has been a growing interest in solving problems in which 
structures interact with fluids. Within the Navy, for example, application 
areas include underwater vibrations, ship silencing, shock response of ships, 
underside slamming, flow- induced vibrations, and the motions of liquids and 
gases in containers. The last two problem areas are also of interest to the 
aerospace community. 

Concern here is with the dynamic structural response of submerged 
structures, including the determination of natural frequencies and the 
response of submerged structures to sinusoidal and general transient 
excitations. In the latter category, a very important application is the 
shock response of submarines to underwater explosions. 

In all problems considered here, the structure is assumed to be initially 
at rest with respect to the fluid; i.e., there is no fluid flow initially. 
Also, the fluid is assumed to behave as an acoustic medium in which each 
material point undergoes only infinitesimal displacement from the static 
equilibrium position. 

Such problems were apparently first formulated by using finite elements 
by Zienkiewicz and Newton (refs. 1 and 2). Choosing pressure as the basic 
fluid unknown, they applied a variational process to obtain the finite-element 
equations for the fluid. These were then coupled to the usual structural 
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equations in a consistent manner. The application of Zienkiewicz's approach 
with NASTRAN (refs. 3 and 4) is straightforward but requires program modifi- 
cations. In fact, a similar approach was adopted in the development of the 
NASTRAN hydroelastic analyzer (refs. 4 and 5). This latter capability is J 
restricted to the case of fluids contained in axisymmetric containers, although 
the motions of the coupled fluid-structure system need not be axisymmetric. 


Another approach to the interaction problem uses fluid displacements 
(rather than pressures) as the fundamental unknowns (ref. 6). In brief, this 
approach converts standard NASTRAN elements into "mock” fluid elements by 
defining (in three dimensions) the 6x6 material matrix G to be 
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In equation (1) , the bulk modulus k is given by 

. 2 
k - pc 


( 2 ) 


where p and c are the unperturbed mass density and sonic speed, respectively, 
of the fluid. To the NASTRAN user, Kalinowski f s displacement approach 
(ref, 6) has the advantages of using a standard version of the program and of 
generating symmetric matrices (in contrast to Zienkiewicz’s approach which 
generates nonsymmetric matrices). However, Kalinowski’s method has the dis- 
advantage of requiring three unknowns per fluid point rather than one 
(pressure) • This penalty affects adversely both the order and the bandwidth of 
the assembled matrices. 


The purpose of the remainder of this paper is to describe two additional 
methods for solving fluid-structure interaction problems. Both use only 
standard versions of NASTRAN. 


A PRESSURE ANALOG METHOD 


All approaches described in this paper treat the fluid as an acoustic 
medium: a compressible, inviscid fluid which is initially stationary and 

undergoes only small amplitude motion, and whose pressure satisfies the wave 
equation 

V 2 p = p/c 2 (3) 

where overdots denote partial differentiation with respect to the time t. The 
boundary condition at the fluid-structure interface can be obtained from 
momentum and continuity considerations: 


k20 


8p 

| 3n = " pU n < 4 > 

rhere n is the unit outward normal from the solid at the fluid-solid inter- 
face, and p is the fluid mass density. A special case of equation (4) occurs 
ft rigid walls where 

| U = 0 < 5 > 

a free surface, in the absence of surface waves, the boundary condition is 
Simply 

P - 0 (6) 

I 

I The goals here are to solve fluid-structure interaction problems in 
i*hich the fluid is described by equations (3) to (6) and to use only stand- 
ard capabilities available in structural analysis codes such as NASTRAN. 

In classical elasticity theory (ref. 7), the x component of the momentum 
aquation is 

I a -t*a +a +pf= pu (7) 

I XX, X xy,y xz,z y x K v J 

where p is the mass density of the (solid) medium, u is the x component of 
displacement, f x is the x component of body force, and the partial 
differentiation with respect to a coordinate is abbreviated with commas: 

< >>x E |;< > < 8 > 

In programs such as NASTRAN, the stress and strain components are generally 
described with vector notation as follows : 


u, + v, 

y x 

v, + w, 

z ’y 

w, + u, 

x z 


where u, v, and w are the three Cartesian displacement components. Thus, for 
linear materials, Hooke's law can be written in terms of a systematic 6x6 
material matrix G, where 

I 

a = G e (101 


By substituting the strain definitions (9) into (10) and (10) into (7) , 
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the wave equation (3) can be obtained if 



u = 

P 

( 11 ) 


V = w 

= 0 

( 12 ) 
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f = 

X 

0 

( 14 ) 


In other words, the x component of displacement u can represent fluid 
pressure p if v and w are fixed everywhere in the fluid, if no body forces 
are applied, and if six of the material constants are as prescribed in equa- 
tion (13). Although the other 15 constants can be arbitrarily chosen, it is 
convenient to choose them so that G is invariant under a coordinate system 
rotation (i.e., G is isotropic). A necessary and sufficient condition for (3 
to be isotropic is that it have the general form 
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(15) 


where X and y are the Lame constants (ref. 7). The only (3 matrix satisfying 
both equations (15) and (13) is 


G = 


2 

pc 


1 
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(16) 


With (3 isotropic, there is no need to be concerned about using finite elements 
whose material matrices are based on some local element coordinate system. 


To summarize, the fluid is modeled with standard elastic finite elements 
(e.g., QDMEM in 2-D and general solids in 3-D) having material properties give 
by equation (16) , where fluid pressure is represented by the x displacement u. 
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j Boundary Conditions 

[interface Condition 

The boundary condition that the pressure p must satisfy at a fluid-solid 
interface is given by equation (4). To evaluate the left-hand side of equa- 
tion (4) , compute the directional derivative of p in the direction of the unit 
outward normal v_ from the fluid at a surface point. Replacing p with its 
.structural analog u yields 


3u n 

— = Vu* v = 
3v — — 


u, v + u, 

X X 


y y 


+ u. 


z z 


(17) 


I By using the constitutive equations (10) and (16) with equation (17) yields 


i 
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(18) 


! 

where the p 
[vector t(v) 


arenthetical expression is equal to the x component of the stress 
acting on a surface with unit outward normal (ref. 7). Hence, 


T (v) 

3u _ x 
3v 2 


(19) 


pc 


If the surface is discretized by a finite number of grid points, the surface 
I traction T( v ) can be replaced by its lumped equivalent 


T (v) _ !x 
x A 


( 20 ) 


where F x is the x component of the force applied to a particular point (on the 
surface with outward unit normal v) to which the area A has been assigned. 
Hence, from equations (19) and (20), the final expression for the directional 
derivative is obtained as 


3u 

3v 


F 

x 


pc A 


( 21 ) 


Since all symbols in equation (21) refer to the fluid, including u which 
represents the pressure p, merely combine equation (21) with the boundary con- 
dition of equation (4) at the interface to obtain the lumped interface condition 


F P = 
x 


(pc) Au 


n 


( 22 ) 


! in which n = - v. and the superscript p has been placed on F x to emphasize that 
the "force" is applied to the pressure variable. That is, if the outward normal 
component of structural acceleration at a point is ii n , the effect on the fluid 
pressure is that of a "force" given by equation (22) applied to the fluid 
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variable p. In. equation (22), p and c refer to the fluid. 

Equation (22) provides the influence of the structural motion on the 
pressure. The inverse relationship (that of fluid pressure on the structural 
motion) is obtained merely by applying a normal force to the structural point 
equal to pA. 


The final set of matrix equations takes the form 


MU + KU « F 


(23) 


where, with the subscripts s and f denoting solid and fluid. 


U 



/^s E J4s 4fi 
-4 S 4fJ — 


o K, f . 


(24) 


In equation (23) , the vector of unknowns U includes both the structural dis- 
placements and the fluid pressures. In equation (24), M fg is a matrix which I 
is assembled by placing the term -(pc) 2 A in the row corresponding to each inter 
face p variable and the column corresponding to the associated structural normaJ 
displacement. Similarly, Kg f is a matrix which is assembled by placing +A in ; 
the row corresponding to each surface normal displacement and the column corre- 
sponding to the associated p. In NASTRAN, both M^ s and Kf s can be inserted by 
using direct matrix input (DMIG cards). 


Infinite Fluids 

The foregoing derivations apply directly to a wide variety of fluid-solid 
systems of finite extent. For structures in "infinite" fluids, interface 
disturbances travel far without reflection so that radiation damping may be 
significant. Although one should generally consider other approaches to solve 
such problems, approximations involving the current method can also be 
derived. 

The approach here is to truncate the fluid model at a distance 
"sufficiently far" from the structure and apply the radiation condition 
(ref. 8) there: 


lE. = _ I iLR (25) 

8n c dt 

This condition has also been referred to as the nonreflecting boundary con- 
dition (ref. 1) or the wave-absorbing boundary condition (refs. 6 and 9). 

Equation (25) can be applied with NASTRAN by combining with equation (21) 3 
where again the pressure p is replaced by its structural analog u. Hence, 
the radiation condition becomes 


T?P * * 

F x = - p cAp (26) 

In other words, a grounded scalar dashpot with damping constant pcA should be 
connected to each p variable. 


Numerical Example 


These concepts can be illustrated by considering the two-dimensional 
problem of a steel ring vibrating in water (ref. 10). The ring has a radius 
of ten inches and a radial thickness of one inch. The inside is evacuated. 
The fluid region is a circular annulus with an outer radius of 32.6 inches. 
The outer surface is a "free" surface (where p=0) . The problem was solved* 
(ref. 10) by modeling with NASTRAN a 90° sector comprised of eight BAR 
elements representing the ring and an 8x6 mesh of QDMEM elements representing 
the fluid. Typical results for the in-fluid natural frequencies (in radians/ 
second) are as follows: 


Fourier NASTRAN Analytical 

Harmonic frequency frequ ency 

0 4333. 4364. 

2 231.6 226.1 

4 1486. 1432. 

6 3884. 3724. 


The analytical results were derived by Schroeder and Marcus (ref. 11) based on 
a method developed by Junger (ref. 12). Calculations reported in Ref. 10 

indicate that grid refinement does result in convergence to the analytical 
results. 


DECOUPLING APPROXIMATIONS 


| For problems in which a structure interacts with a fluid of infinite 
f extent, it may be expensive to model enough fluid to produce satisfactory 
results, even if the radiation condition (eqs. (25) and (26)) is applied at 
some arbitrary outer surface. Since the interest is in the structure rather 
(than in the fluid, various schemes can also be used that approximate the fluid 
effects without actually modeling the fluid (refs. 13 and 14). These schemes, 
generally referred to as decoupling approximations, result in an analytical 
expression (in differential equation form) describing the relationship between 
the fluid pressure at the interface and the interface motion. 

For certain dynamics problems in which the source of the disturbance is 
in the fluid rather than in the structure (e.g., underwater explosions), the 
decoupling approximation supplies only the radiation pressure. However, the 
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other contributors to the total dynamic pressure (the incident free field 
pressure and the scattered pressure) can be computed as if the structure were 
rigid and stationary and thus comprise the usual right-hand-side forcing 
function. 


One decoupling approximation which has been used with NASTRAN was 
developed (but not yet published) by Dr. Hansen Huang of the Naval Research 
Laboratory. For axisymmetric cylinders modeled with conical shell elements, 
the interface pressure-motion equation is approximated by 


p n + (cn/r)p n = pcw n 


(27) 


where p 


c 

P 

r 


w 


n 


is the n th Fourier harmonic of the radiation pressure 
is the speed of sound in fluid 
is the mass density of fluid 
is the radius of cylinder 

is the n th Fourier harmonic of the outward normal component of 
shell displacement 


Equation (27) is equivalent to the so-called doubly asymptotic approximation 
of Geers (ref. 13), which was formulated for more general three-dimensional 
situations . 


To illustrate how approximations such as equation (27) are applied with 
NASTRAN consider a cylindrical shell subjected to an underwater shock loading 
( re f t 15). In the absence of radiation pressure p, the usual transient matrix 

equation would apply 


mu + ku = f_(t) 


(28) 


where the forcing function f consists of the sum of the incident and scattered 
pressures. The radiation pressure, which depends on the shell motion, 
supplies an additional load on the structure so that the equilibrium equations 
for the structure, from equation (28), become 


m u 


u + k u + h£ = f_(t) 


(29) 


where £ is the vector of radiation pressures at each point, and h is a matrix 
of area factors converting pressure to force. 


Because of the mutual dependence of the shell motion and the radiation 
pressure on each other, equations (27) and (29) must be solved simultaneously. 
This solution is accomplished in NASTRAN by defining a new set of scalar degree 
of freedom representing the radiation pressure at each surface point. The set 
of all such pressures is grouped into the vector jd. Thus, to NASTRAN, the 
vector of unknowns (displacements and pressures) is 


k26 


(30) 



After equation (27) has been written at each point and combined with equa- 
ion (29) , the resulting matrix equation which NASTRAN integrates is of the 
isual form 

I 

MU + B U + K U = F(t) (31) 

?he contributions to M, JB, and K corresponding to the u variables are computed 
lirectly by NASTRAN based on geometry and material properties. The contribu- 
tions to M, B, and K corresponding to the radiation pressure variables (p) are 
.he coefficients appearing in equation (27). Hence, they are computed exter- 
nally to NASTRAN and input directly (using the TF bulk data card) into the 
Program. Thus, the matrices M, B_, K, and F in equation (31) take the unsym- 
jietric form 


“ m 0“ 


1 

o 

o 

1 

-fit 

( i (t) ) 



B = 


K = 

F - { > 

(32) 

_d 0_ 


Lo a J 

Lo bj 

\ 0 f 



NASTRAN integrates equations (32) by using the Newmark-Beta finite-difference 
algorithm, which is unconditionally stable if equation (32) is stable. 

Unfortunately, specific numerical results are not yet available for 
general publication. However, early success with decoupling approximations 
nas convinced the authors that such approaches have great promise as an alter- 
native to the explicit finite element modeling of fluids of "infinite" extent. 


GENERAL REMARKS 

I 


Two methods by which NASTRAN can be applied without modification to the 
solution of fluid-structure interaction problems have been described. 

The first, the structural analog, is the "lumped" equivalent of the 
consistent formulation described by Zienkiewicz and Newton (refs. 1 an d 2). 
lowever , the lumped approach has the advantage for the NASTRAN user that it 
pan be. applied without program modification. Otherwise, the two pressure 
formulations have similar characteristics, including nonsymmetric matrices. 
This is in sharp contrast to displacement approaches (ref. 6) which assemble 
symmetric matrices, but at the expense of having three times the n umb er of 
fluid unknowns. 

It should be emphasized that although the analog method was developed 
specifically for an acoustic fluid, the same technique could be applied to a 
pri.de variety of problems in mathematical physics, including heat conduction 
(ref. 16), the Helmholtz equation, electrical or magnetic potential problems, 
the torsion of prismatic bars, potential fluid flow, or seepage through porous 
media. In the present context of fluid-structure interaction problems, the 
calculation of added mass matrices can be accomplished directly, since it 
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involves solving Poisson's equation in the fluid region where the source 
terms occur only at the fluid-structure interface. Such equations can also bd 
solved with standard three-dimensional steady-state heat conduction codes 
such as those contained in NASTRAN Level 15.5 (ref. 4), the Navy's thermo- i 
structural NASTRAN (refs, 17 and 18), or the CINDA thermal analysis program 

(ref. 19). 

For infinite media, the explicit modeling of the fluid is often 
uneconomical compared with some of the competitive methods. For fluid- 
structure interaction problems, the use of decoupling approximations pro- 
vides an attractive alternative. For other general field problems (e.g., 
potential flow), integral equation techniques are widely used (ref. 20). 
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